%% process_plot_scalar
function process_plot_scalar(mesh, fespace, SOL, varargin)
x = fespace.node(:,1);
y = fespace.node(:,2);
if isempty(varargin)
    SOL = full(SOL);
    figure("WindowStyle", "docked");
    elem = convert_elem(mesh, fespace);
    trisurf(elem, x, y, SOL);
    colorbar;
    xlabel("x");
    ylabel("y");
    zlabel("uh");
    axis("equal", "tight");
    view(-45,60);
else
    figure("WindowStyle", "docked");
    scatter3(x, y, SOL, 20, SOL, "filled");
    colorbar;
    xlabel("x");
    ylabel("y");
    zlabel("uh");
    axis("equal", "tight");
    view(-45,60);
end
end

%% convert_elem
function new_elem = convert_elem(mesh, fespace)
if strcmp(fespace.basis_type, "P1")
    new_elem = fespace.elem;
elseif strcmp(fespace.basis_type, "P1b")
    mod = [1, 2, 4;
        2, 3, 4;
        3, 1, 4];
    new_elem = reshape(fespace.elem(:, mod(:)), 3*mesh.Ne, 3);
elseif strcmp(fespace.basis_type, "P2")
    mod = [1, 4, 6;
        2, 5, 4;
        3, 6, 5;
        4, 5, 6];
    new_elem = reshape(fespace.elem(:, mod(:)), 4*mesh.Ne, 3);
elseif strcmp(fespace.basis_type, "P3")
    mod = [1, 4, 9;
        2, 5, 7;
        3, 6, 8;
        4, 7, 10;
        5, 8, 10;
        6, 9, 10;
        7, 5, 10;
        8, 6, 10;
        9, 4, 10];
    new_elem = reshape(fespace.elem(:, mod(:)), 9*mesh.Ne, 3);
elseif strcmp(fespace.basis_type, "P4")
    mod = [1, 4, 12;
        2, 5, 10;
        3, 6, 11;
        4, 7, 13;
        5, 8, 14;
        6, 9, 15;
        7, 10, 14;
        8, 11, 15;
        9, 12, 13;
        10, 5, 14;
        11, 6, 15;
        12, 4, 13;
        13, 7, 14;
        14, 8, 15;
        15, 9, 13;
        13, 14, 15];
    new_elem = reshape(fespace.elem(:, mod(:)), 16*mesh.Ne, 3);
elseif strcmp(fespace.basis_type, "P5")
    mod = [1, 4, 15;
        2, 5, 13;
        3, 6, 14;
        4, 7, 16;
        5, 8, 17;
        6, 9, 18;
        7, 10, 19;
        8, 11, 20;
        9, 12, 21;
        10, 13, 17;
        11, 14, 18;
        12, 15, 16;
        13, 5, 17;
        14, 6, 18;
        15, 4, 16;
        16, 19, 21;
        16, 7, 19;
        17, 8, 20;
        18, 9, 21;
        19, 10, 17;
        20, 11, 18;
        21, 12, 16;
        17, 20, 19;
        18, 21, 20;
        19, 20, 21];
    new_elem = reshape(fespace.elem(:, mod(:)), 25*mesh.Ne, 3);
end
end

